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I. 


SCOPE  OF  PROGRAMS 


As  indicated  in  References  2  and  3,  the  analyses  by  Dr,  Frank  Lane 
proved  to  be  acceptable,  using  45  degrees  of  freedom,  only  for  flat  plates, 
twisted  bars  ,  and,  in  general,  for  any  geometrical  configuration  which 
was  not  characterized  by  severe  twist,  and  large  overhang  of  thin  leading 
and  trailing  edges.  Since  many  propeller  blades  do  exhibit  the  geometrical 
configuration  described  above,  it  was  suggested  that  the  number  of  degrees 
of  freedom  be  increased,  to  accommodate  all  propeller  blade  configurations. 

Since  the  computer  machine  time  increases  roughly  as  the  square  of 
the  number  of  degrees  of  freedom,  it  was  decided  to  re-code  the  programs 
described  in  References  2  and  3,  for  the  IBM  7090  computer,  instead  of 
for  the  IBM  704.  Due  to  the  size  of  the  system,  it  was  necessary  to  sub¬ 
divide  the  program  into  4  CHAIN  LINKS,  to  accommodate  the  storage 
requirements.  During  the  course  of  field  testing,  it  was  found  that  single¬ 
precision  arithmetic  was  inadequate  in  one  section,  and  therefore  it  was 
necessary  to  utilize  triple -precision  arithmetic. 

To  illustrate  the  sharp  contrast  in  scope,  between  the  subject  programs 
and  their  IBM  704  predecessors,  up  to  108  degrees  of  freedom  may  now  be 
utilized,  as  compared  with  45;  10  tapes  versus  3,  up  to  146  mesh  points  versus 
56  (integrating  quadrature),  and  about  1  hour  of  IBM7090  computer  time 
versus  40  minutes  of  IBM  704  time.  As  before,  provision  is  made  to  restart 
the  program  with  a  minimum  loss  of  machine  time,  in  the  event  it  must  be 
removed  from  the  computer  prematurely  due  to  lack  of  available  time.  In  ^ 
addition,  the  program  may  be  easily  restored  in  the  event  of  a  machine  error. 


2. 


No  attempt  will  be  made  in  this  report  to  repeat  Dr.  Lane's  analyses, 
which  are  described  in  Reference  1,  and  summarized  in  References  2  and  3. 
Only  that  information  which  pertains  to  these  IBM  7090  programs  will  be 
included.  The  expanded  sets  of  Rayleigh  functions  are  listed  in  Appendix  I. 
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PROPELLER  BLADE  VIBRATION  ANALYSIS 


The  purpose  of  this  IBM  computer  program  is  to  generate  the  natural 
frequencies,  mode  shapes,  and  modal  stress  distributions  for  the  ten  lowest 
modes  of  vibration  of  propeller  blades.  This  information,  coupled  with  the 

fact  that  the  eigenvectors  are  normalized  with  respect  to  the  strain  energy 
T 

matrix,  (q  Vq  *1),  is  sufficient  to  conduct  any  response  analyses  of  propeller 
blades. 

After  accepting  the  input  data  which  describe  the  blade  geometry,  the 
program  generates  a  polynomial  which  fits,  in  the  least  squares  sense,  the 
camber  distribution  over  the  planform  of  the  blade'.  The  program  then 
constructs  a  quadrature  mesh  which  is  later  used  to  integrate  the  strain  energy 
and  kinetic  energy,  matrices  over  the  blade  planform.  The  two-matrix  eigen¬ 
value  problem  is  reduced  to  a  one  matrix  eigenvalue  problem  by  factoring  the 
strain  ehergy  matrix  in  triple-precision.  This  problem  is  then  solved,  and  the 
resulting  eigenvalues  and  eigenvectors  are  written  as  output.  The  relative 
stresses  and  displacements  are  calculated  at  mesh  points  over  the  blade  plan- 
form,  and  are  written  as  output. 

This  program  will  process  sets  of  input  data  automatically.  Core  memory 
is  salvaged  on  tape  three  times  during  the  execution  of  the  program  which  permits 
the  program  to  be  restarted  with  a  minimum  loss  of  time  in  the  event  of  a 


machine  failure. 
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This  program  has  been  written  for  a  32K  (32,  768  words)  capacity 
IBM  7090  computer,  and  utilizes  ten  tape  units  , 

CHAIN  LINK  1 

The  first  two  input  cards,  identification  card  and  the  first  data  card, 
are  read  and  immediately  written  on  tape  A3.  The  program  then  reads 
and  writes  the  remainder  of  the  input  data  and  calculates  the  thickness  and 
the  local  distance  from  the  y*,  y^  plane  of  each  input  point.  The  distribution 
of  camber,  of  the  propeller  blade,  is  then  approximated  with  a  least-squares 
surface  fit.  The  LeGendre-Gauss  Quadrature  weighting  constants  are  stored 
and  the  coordinates  of  the  quadrature  mesh  points  are  calculated.  The 
program  then  dumps  core  memory  on  tape  A5  and  calls  for  CHAIN  LINK  2. 
CHAIN  LINK  2 

This  section  generates  the  potential  (strain)  energy,  and  the  kinetic  energy 
matrices.  At  each  quadrature  mesh  point,  the  program  generates  the  funda¬ 
mental  forms,  the  Rayleigh  functions,  and  the  elements  of  the  T...  and  N 

(i) 

matrices,  as  well  as  the  local  thickness.  Utilizing, these  values,  the  local 
contribution  to  the  strain  energy  and  kinetic  energy  matrices  is  formed  and 
multiplied  by  the  local  quadrature  weight.  All  underflows  are  set  to  zero. 

The  program  scans  the  quadrature  mesh  points,  adding  the  weighted  contri¬ 
butions.  The  integrals  are  then  modified  to  account  for  the  root  curvature. 
After  the  integrals  are  formed,  core  memory  is  again  salvaged  on  tape  AS, 


and  CHAIN  LINK  3  is  called. 
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CHAIN  LINK  3 

The  two-matrix  eigenvalue  problem  is  expressed  as  Tq  =  ~  Vq,  Using 

A 

T 

triple -precision  arithmetic  the  strain  energy  matrix  is  factored,  V=M  M, 

-T  i. 

where  M  is  an  upper  triangular  matrix.  Then  the  system,  M  TM 
1  -T  - 1 

or  Fy  ■  ~.y  ,  where  F=M  TM  and  y  =  Mq,  is  formed.  This  eigenvalue 

A 

problem  is  solved  using  the  SHARE  routine,  NYEVV,  and  the  true  eigenvectors 
calculated,  q  =  M  V*  These  eigenvectors  are  normalized  with  respect  to  the 
strain  energy  matrix  (q  Vq  ■  I).  After  saving  core  memory  on  tape  A5,  the 


program  calls  for  CHAIN  LINK  4. 


CHAIN  LINK  4 

The  ten  lowest  eigenvalues  and  their  corresponding  eigenvectors  are  stored 
in  ascending  order,  and  the  mesh  points  at  which  the  displacements  and  stresses 
will  be  evaluated,  are  calculated. 

After  writing  the  eigenvalues  and  eigenvectors  for  the  lowest  ten  modes,  and 
the  mesh  point  coordinates  on  output  tape  A3,  the  program  calculates  the  funda¬ 
mental  forms,  Rayleigh  functions,  elements  of  the  T  and  N,.,  matrices,  the 

(i)  M 

thickness  at  each  mesh  point,  and  then  calculates  displacements  and  stresses. 

Since  the  mesh  is  very  dense,  stresses  are  calculated  at  alternate  mesh  points 
while  displacements  are  calculated  at  all  mesh  points,  for  all  the  ten  modes  of 
vibration.  There  are  a  total  of  111  mesh  points  including  the  tip  of  the  blade  and 
11  points  along  the  blade  root.  The  stresses  and/or  displacements  at  each  mesh 
point  are  written  on  output  tape  A3.  Control  is  returned  to  CHAIN  LINK  1  to  process 
the  next  run,  if  any.  If  no  more  data  is  available,  the  program  exits. 

The  approximate  computer  production  time  for  each  set  of  input  cards  for 
an  8  x  9  quadrature  mesh  and  90  degrees  of  freedome  is  approximately  1  hour. 
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PROPELLER  BLADE  STATIC  STRESS  ANALYSIS 
A.  Program  for  Distributed  Loads 

The  purpose  of  this  IBM  computer  program  is  to  generate  the 
displacement  and  stress  distribution  of  propeller  blades  for  prescribed 
pressure-jump  loadings.  Provision  is  made  for  a  "peaked"  distributed  loading 
with  the  singularity  at  the  leading  edge.  The  stresses  are  given  directly; 
the  displacements  given  must  be  divided  by  Young's  Modulus. 

After  accepting  the  input  data  which  describes  the  blade  geometry,  tlve 
program  processes  this  data,  and  generates  functions  which  accurately 
represent  the  blade  configuration.  A  quadrature  mesh  is  constructed  which 
is  then  used  to  integrate  the  strain  energy  matrix  over  the  blade  planform. 

The  program  reads  the  pressure  jump  load  data,  and  integrates  this  informa¬ 
tion  over  the  blade  planform  to  generate  the  forcing  vector.  To  solve  the 
resulting  system  of  linear  equations,  the  strain  energy  matrix  is  factored  in 
triple-precision,  and  its  inverse  is  computed.  This  approximate  inverse  is 
improved  (if  possible),  and  the  solution  vector  is  calculated  and  improved  by 
use  of  iteration  technique?.  See  Appendix  II.  The  stresses  and  displacements 
at  mesh  points  over  the  blade  planform,.  are  then  calculated. 

This  program  will  process  any  sequence  of  inputs  consisting  of  geometrical 
input  data  describing  a  particular  blade,  followed  by  any  number  of  loading 
conditions  for  that  blade,  followed  by  another  blade,  etc.  All  inputs  are 
written  on  the  output  tape,  so  that  any  error  maybe  detected,  and  all 
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output  listings  uniquely  identified.  Core  memory  is  salvaged  on  tape  four  times 
during  the  execution  of  the  program  which  permits  the  program  to  be  restarted 
with  a  minimum  loss  of  time  in  the  event  of  a  machine  failure. 

This  program  has  been  written  for  a  32K  (32,  768  words)  capacity  IBM 
7090  computer,  and  utilizes  ten  tape  units. 

The  program  is  divided  into  four  CHAIN  LINKS,  See  pages  32  to  38 
for  detailed  flow  diagrams. 

CHAIN  LINK  1 

The  first  two  input  cards,  identification  card  and  the  first  data  card, 
are  read  and  immediately  written  on  tape  A3.  The  program  then  reads 
and  writes  the  remainder  of  the  input  data  and  calculates  the  thickness  and 
the  local  distance  from  the  y*,  y^  plane  of  each  input  point.  The  distribution 
of  camber  of  the  propeller  blade,  is  then  approximated  with  a  least-squares 
surface  fit.  The  LeGendre -Gauss  Quadrature  weighting  constants  are  stored 
and  the  coordinates  of  the  quadrature  mesh  points  are  calculated.  The 
program  then  dumps  core  memory  on  tape  A5  and  calls  for  CHAIN  LINK  2. 
CHAIN  LINK  2 

This  section  generates  the  potential  (strain)  energy  matrix.  At  each 
quadrature  mesh  point,  the  program  generates  the  fundamental  forms,  the 
Rayleigh  functions,  and  the  elements  of  the  T^^j  and  N^^^^  matrices,  as  well 
as  the  local  thickness.  Utilizing  these  values,  the  local  contribution  to  the 
strain  energy  matrix  is  formed  and  multiplied  by  the  local  quadrature  weight. 
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All  underflows  are  set  to  zero.  The  program  scans  the  quadrature  mesh 
points,  adding  the  weighted  contributions.  The  integral  is  then  modified  to 
account  for  the  root  curvature.  After  the  integral  is  formed,  core  memory 
is  again  saved  on  tape  A5,  and  CHAIN  LINK  3  is  called. 


CHAIN  LINK  3 

A  code  word  is  read  and  written  which  describes  the  type  of 

blade  in  the  event  the  loading  is  "peaked"  at  the  leading  edge.  The  pressure- 

jump  distributed  loading  is  then  read  and  written.  A  forcing  vector  is  formed 

and  integrated  over  the  blade  planform. 

Core  memory  is  saved  on  tape  A5.  The  strain  energy  matrix  is  then 

T 

factored  into  triangular  matrices  in  the  form,  V  =  M  M,  using  triple-precision 

arithmetic,  and  the  matrix  M  ^  is  generated.  The  inverse  is  then  simply, 

-1  -1  -T 

V  =  M  M  ,  and  is  written  on  binary  tape,  B2,  for  subsequent  loading 
conditions,  if  desired.  The  solution  vector  is  then  calculated  and  improved  by 
iteration.  All  subsequent  loading  conditions  use  this  inverse;  only  the  new  forcing 
vector  is  computed.  After  saving  core  memory  again  on  tape  AS,  the  program 


calls  CHAIN  LINK  4. 
CHAIN  LINK  4 


The  mesh  points  at  which  the  displacements  and  stresses  will  be 
evaluated,  are  calculated.  After  writing  the  solution  vector  and  the  mesh 
points  on  the  output  tape,  the  program  calculates  the  fundamental  forms, 
Rayleigh  functions,  elements  of  the  T^.^  and  N^.^  matrices,  the  thickness, 
and  then  calculates  displacements  and  stresses  at  each  mesh  point.  These 
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values  are  written  on  the  output  tape.  There  are  a  total  of  111  mesh  points 
including  the  tip  of  the  blade  and  11  points  along  the  blade  root.  If  a  new 
loading  condition  is  to  be  read,  CHAIN  LINK  3  is  called.  If  a  new  blade  is  to 
be  analyzed,  control  is  transferred  to  CHAIN  1(INK  1.  Otherwise,  the  program 
exits . 

The  approximate  computer  production  time  for  each  complete  set  of 
input  data  cards  for  an  8  x  9  quadrature  mesh  and  90  degrees  of  freedom 
is  approximately  1  houri  .  For  each  subsequent  set  of  pressure-jump  inputs 
alone,  the  computer  production  time  is  approximately  1/2  hour. 
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n  GEOMETRY 

Figure  1  illustrates  the  propeller  blade  geometry  and  the  coordinate 
systems  utilized.  The  (y*,y*)  coordinates  form  a  two-dimensional  Cartesian 
system  in  the  transverse  projection  at  the  blade  and  serve  as  a  Gaussian 
coordinate  system  in  the  blade  middle-surface.  The  coordinate  y’  is 
normal  to  the  local  blade  middle-surface,  directed  so  that  (y*,y*,  y*)  form 
a  right-hand  system.  It  is  to  be  noted  that  the  camber  and  twist  of  the 
blade  result  in  non-orthogonality  of  the  (y',y^)  coordinate  system.  The 
third  coordinate  y’  is,  however,  locally  normal  to  y*  and  y*. 

The  blade  middle -surface  is  defined  completely  by  describing  its 
projected  planform  in  the  (y*,y*)  plane  and  its  local  distance  f(y',y*)  from 
this  plane  as  shown  in  Figure  1.  The  description  of  the  thickness  distribution 
serves  to  complete  the  geometric  description  of  the  blade. 

The  thickness  at  any  point  is  found  by  linear  interpolation  qf  the  input 

data. 
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in  INPUT  FORMAT  -  VIBRATlfi|N 


Card 

1 

2 

3.6,9,  ... 
4,7,  10,... 
5,  8,  1 1,. .  . 
Final 


Input  Information 
Identification 

Poisson's  Ratio,  Blade  Root  Geometry,  Siae  of  System 

Section  Properties 

Offsets  to  Back  Face  Along  Section 

Offsets  to  Front  Face  Along  Section 

Properties  at  Blade  Tip 


Refer  to  Figure  2. 

Card  No.  1: 

Any  identification  statement  of  72  characters  or  less  will  be  accepted. 
The  statement  may  be  punched  in  columns  2  -  72  of  a  card;  the  digit,  "1" 
should  be  punched  in  column  1  for  output  convenience. 

Card  No.  2: 


Data 

Columns 

Format 

V 

1-5 

x.xxx 

^R 

7-12 

XX.  XXX 

14-19 

XX,  XXX 

R 

21-26 

XX.  XXX 

R 

28-33 

XX.  XXX 

N 

37-38 

XX 

V 
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The  first  word,  V,  is  Poisaon*'8  Ratio.  The  next  4  words,  ^R’  ^li 

are  indicated  in  Figure  2.  The  number  of  cjegrees  of  freedom  for  each 
displacement  component  is  N.  Thus  for  63  degrees  of  freedom,  N  is  21. 

If  a  full-scale  problem  is  being  run  (108  degrees  of  freedom),  N  may  be  left 
blank. 

Card  No.  3; 


Data 

Columns 

Format 

%  Radius 

1-2 

XX 

Radius 

4-9 

XX.  XXX 

or  XXX.  XX 

Span 

11-16 

XX.  XXX 

or  XXX.  XX 

Pitch 

18.24 

+XX.XXX  or  +XXX.XX 

Rake 

26-32 

iX.XXXX 

Skew 

34-40 

+X.XXXX 

Card  No.  4: 

Data\ 

Colurms 

Format 

A 

1-6 

+X.XXX  or  +XX.XX 

E 

7-12 

F 

13-18 

i 

1 

G 

19-24 

1 

1 

H 

25-30 

J 

31-36 

K 

37-42 

L 

43-48 

M 

49-54 

N 

55-60 

P 

61-66 
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Card  No.  5: 

Identical  in  format  to  Card  No.  4. 
Final  Card 


Data 

Columns 

Format 

- 

1-2 

Must  be  blank 

Radius 

4-9 

XX.  XXX  or  XXX.  XX 

Thickness 

11-16 

XX.  XXX 

Pitch 

18-24 

XXX.  XXX 

Rake 

26-32 

+x.  xxxx 

m 

Skew 

34-40 

+X.XXXX 

I 

H  The  program  can  accept  a  maximum  of  12  sections  corresponding  to  a 

total  of  40  input  cards  per  set. 

I  Cards  4  and  5  specify  the  thickness  offsets  along  a  section  at  equal  chord 

!  increments,  from  leading-edge  to  trailing-edge.  The  offset  of  A  is  the 

offset  near  the  leading  edge,  the  offset  of  E  is  at  one-tenth  the  chord 
distance  from  the  leading  edge,  etc.  The  offset  of  P  refers  to  the  trailing- 
edge.  See  sketch  below. 


INPUT  FORMAT  -  STRESS. 


Card 


1 

2 

3.6.9.. .. 
4,  7,  10,. . . 

5.8.11.. .. 


T-l 


T 

T+1,  T+3,  T+5, .  . . 


T+2,T+4,  T+6,  V  .. 


C 


Input  Information 
Identification 

Poisson's  Ratio,  Blade  Root  Geometry 

Section  Properties 

Offsets  to  Back  Face  Along  Section 

Offsets  to  Front  Face  Along  Section 

Properties  of  Blade  Tip 

Code  Card 

Pressure-jumps  at  points  on  sections 
corresponding  respectively  to  those  sections 
of  input  cards  3,  6,  9,  . . .  These  points  are 
from  leading  edge  to  25%  of  chord. 

Pres  sure -jumps  at  points  on  sections 
corresponding  respectively  to  those  sections 
of  input  cards  3,  6,  9,  .  .  .  These  points  are 
from  30%  of  chord  to  trailing  edge. 

Control  Card 


Next  set  of  inputs,  if  any. 

Refer  to  Figure  2. 

Card  No.  1: 

Any  identification  statement  of  71  characters  or  less  will  be  accepted. 

The  statement  must  be  punched  in  columns  2-72  of  a  card.  The  digit,  "1" 
i-i 

should  be  punched  in  column  1  for  convenience  of  output. 
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Card  No.  2; 


Data 

Columns 

Format 

U 

1-5 

x.xxx 

w 

7-12 

XX.  XXX 

R 

14-19 

XX.  XXX 

21-26 

XX.  XXX 

28-33 

XX.  XXX 

N 

37-38 

XX 

The  first  word,  V,  is  Poisson's  Ratio.  The  next  4  words, 

are  indicated 

in  Figure  2.  The  number  of  degrees  of  freedom  for  each 

displacement  component  is  N. 

Thus  for  63  degrees  of  freedom,  N  is  21. 

If  a  full-scale  problem  is  being 

run  (108  degrees  of  freedom),  N  maybe  left 

blank. 

Card  No.  3: 

i 

Data 

Columns 

F  ormat 

%  Radius 

1-2 

XX 

Radius 

4-9 

XX.  XXX  or  XXX.  XX 

Span 

11-16 

XX.  XXX  or  XXX.  XX 

Pitch 

18-24 

+XX.  XXX  or  ±XXX.  XX 

Rake 

26-32 

+X.XXXX 

Skew 

34-40 

+X.  XXXX 

Card  No.  4: 
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Data 

Columns 

F  ormat 

A 

1-6 

tX.XXX  or  +  XX.  XX 

E 

7-12 

F 

13-18 

G 

19-24 

H 

25-30 

J 

31-36 

K 

37-42 

L 

43-48 

M 

49-54 

N 

55-60 

1 

r 

P 

61-66 

Card  No.  5: 

Identical  in 

format  to  Card  No.  4. 

Card  T-1; 

Data 

Columns 

F  ormat 

- 

1-2 

Must  be  blank 

Radius 

4-9 

XX.  XXX  or  XXX.  XX 

Thickness 

11-16 

XX.  XXX 

Pitch 

18-24 

XXX.  XXX 

Rake 

26-32 

+x.  xxxx 

Skew 

34-40 

+X.XXXX 

Card  T: 

Data 

Column 

Format 

Code 

1 

X 
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Card  T  +  1: 


Data 

•  Columns 

Po 

1-7 

Pi. 25 

9-15 

j  ^2.5 

17-23 

J  ”5.0 

25-31 

1  _ 

P7.5 

33-39 

QQ 

i 

41-47 

1  ^15 

49-55 

■ 

57-63 

1  ”25 

65-71 

1  Card  T  +  2; 

_  Data 

Columns 

1 

P30 

1-7 

1  P40 

9-15 

17-23 

1  50 

1  1 

0^ 

0 

25-31 

1  P70 

33-39 

■1 

OD 

0 

41-47 

■  P 

*^90 

49-55 

1 

57-63 

0 

0 

65-71 

I 

1 


F  ormat 
XXX.  XXX 


Format 
XXX.  XXX 
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Card  C: 

Data  Column  Format 

Code  1  X 

The  program  can  accept  a  maximum  of  12  sections  corresponding 
to  a  total  of  65  input  cards  per  set. 

Cards  4  and  5  specify  the  thickness  offsets  along  a  section  at  equal 
chord  increments,  from  leading-edge  to  trailing-edge.  The  offset  of  A 
is  the  offset  near  the  leading  edge,  the  offset  of  E  is  at  one-tenth 
the  chord  distance  from  the  leading  edge,  etc.  The  offset  of  P  refers  to 
the  trailing  edge.  See  sketch  below. 


Card  T  specifies  the  type  of  peakedloadihg  if  any.  The  code  is  "1" 
for  a  "peaked"  load  and  a  super -cavitating  blade,  0  for  a  "peaked"  load  and 
a  fully  wetted  blade,  and  blank  for  any  load  which  is  not  "peaked". 

Cards  T+l  and  T+2  specify  the  complete  pressure-jump  loading  acting  on 
that  section  described  on  cards  3,  4  and  5.  The  subscripts  of  the  p's  arc  the 
percentages  of  chord,  counting  from  leading  edge  to  trailing  edge.  For  example, 
p^  is  the  pressure-jump  loading  acting  at  the  leading  edge  of  the  section  under 
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consideration;  p  acts  at  20%  of  chord,  or  one-fifth  the  distance  from 
20 

the  leading  edge  (this  corresponds  to  F  on  the  above  sketch).  Of  course, 

PlOo  is  pressure-jump  loading  at  the  trailing  edge  (invariably,  zero). 

Cards  T+3  and  T  -f  4  specify  the  pressure-jump  loading  on  that  section 
described  on  cards  6,  7,  and  8;  cards  T  +  5  and  T  +  6  pertain  to 
cards  9,  10,  and  11,  etc.  For  a  peaked  loading,  P0-P1.25  P2.5 

The  control  card  dictates  the  action  of  the  program  after  a  "run"  has 
been  completed.  If  no  more  inputs  follow,  the  control  card  may  be  blank. 

If  it  is  desired  to  conduct  an  additional  analysis  on  the  same  blade  using  a 
different  loading,  the  control  card  must  have  a  "2"  punched  in  column  1, 
and  be  followed  immediately  by  cards  specifying  the  loading,  as  described 
above  (cards  T  through  C).  Should  it  be  desired  to  conduct  an  analysis  on 
another  blade  (no  matter  what  the  loading  may  be),  a  "1"  must  be  punched 
in  Column  1,  followed  by  cards  describing  the  new  geometry  and  corresponding 
pressure-jump  loading  (cards  1,  2,  3...,C). 

The  formats  indicated  are  representative.  The  decimal  points  may  be 
shifted  in  either  direction,  so  long  as  the  input  "word"  is  confined  within  the 
columns  specified. 

The  program  assumes  that  the  tip  of  the  blade  is  a  point  (span  equal 
to  zero)  as  is  the  case  of  a  propeller  blade.  If  the  section  at  the  blade  tip 
has  some  span,  such  as  the  flat  plates  (Ref.  2  and  3  )  ,  it  will  be  necessary 
to  input  two  sections  extremely  close  to  the  tip.  In  general,  both  sections 
should  be  within  the  outer  1%  of  the  total  blade  span. 
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rV  OPERATING  INSTRUCTIONS 


These  programs  consist  of  rather  large  binary  decks  which  should  be 
loaded  on  tape  using  off-line  peripheral  equipment  such  as  the  IBM  1401,  and 
mounted  as  input  tape  A2  on  the  main  computer.  The  data  cards  are  placed 
in  the  hopper  of  the  on-line  card  reader  as  detailed  below.  . 

The  IB  Monitor  system  must  be  used  to  run  this  program  on  the  IBM 
7090.  The  tape  assignment  cards  (lOU  tables)  are  included  as  part  of  the 
program  deck.  Tapes  A1  through  A5  and  B1  through  B5  are  required  for 
these  programs.  All  output  data  is  written  on  tape  A3  which  is  to  be  printed 
using  program  control.  Tape  A5  is  used  as  a  salvage  tape  to  store  the  entire 
contents  of  core  periodically,  so  as  to  minimize  loss  of  machine  time.  The 
restart  procedure  is  described  below.  Tape  A1  contains  the  IB  Monitor 
Systems  tape,  tape  A4  is  used  to  store  the  chain  links  of  the  program,  and 
tapes  B-1  through  B-5  are  used  as  intermediate  tapes  by  the  Monitor  system 
and  as  binary  scratch  tapes  by  the  program  itself.  All  tapes  should  be  in  the 
high  density  mode  to  conserve  machine  time. 

Restart  Procedure 

A.  In  the  event  the  program  must  be  interrupted,  due  to  some 
scheduling  conflict  or  similar  reason,  dismount  tape  A5,  and  save.  Print  tape 
A3  and  save  this  abridged  output. 

1.  To  restart  the  vibration  program,  re-mount  tape  A5  (but  do  not 
"ready"  the  tape  unit),  the  systems  tape,  Al>,  and  reload  the  program  deck  and  data  in 
the  usual  procedure,  described  above,  as  for  a  new  run.  As  soon  as 
EXECUTION  commences,  (on-line  printer  will  indicate  this  word),  depress 
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MANUAL  key  to  MANUAL  position,  then  "ready"  tapeAS.  Depress  the 
console  keys  to  read,  0  02000  0  00175,  (transfer  to  location  (0175)g).  Press 
ENTER  INSTRUCTION  key.  Depress  MANUAL  key  to  AUTOMATIC  position. 
Press  START  key. 

Tape  A5  will  then  be  read,  and  the  most  recent  condition  of  core, 
prior  to  the  last  salvage,  will  be  restored.  The  program  will  then 
commence  to  execute,  automatically, 

2.  To  restart  the  stress  program.  ,  it  is  necessary  that  the 
pres  sure -jump  data  be  in  core  memory  at  the  time  of  salvage.  Thus,  the 
program  must  have  initially  executed  either  the  3rd  or  4th  dump  on  tape  A5, 
both  taking  place  in  the  third  CHAIN  LINK,  prior  to  manual  termination  of 
run.  If  so,  the  procedure  described  above  for  the  vibration  program,  will 
apply, otherwise,  the  program  must  be  run  from  the  beginning.  Note  that 
sequential  runs  are  not  possible  when  the  program  is  restarted. 

B.  In  the  event  of  a  machine  failure,  the  restart  procedure  is  some¬ 
what  simpler.  Transfer  via  the  console  to  location  (0175)g  as  described 
above,  and  the  program  will  be  restored.  In  this  case,  the  restrictions  for 
the  stress  programs  do  not  apply. 

The  only  STOP  anticipated,  is  when  the  accumulated  round-off  error  causes 
a  number,  whose  square-root  is  required,  to  become  negative.  This  will  only 
occur  in  LINK  3,  and  the  program  will  exit  automatically,  with  a  statement 
describing  the  cause. 

This  situation  can  only  be  remedied  by  re-running  the  program  using  lower 
number  of  degrees  of  freedom.  For  a  full  discussion  of  this  difficulty,  refer 
to  section  of  conclusions. 
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Another  STOP  which  may  occur,  is  at  location  (210)  .  This  is  the 

o 

result  of  being  unable  to  write  on  tape  A5  due  (probably)  to  a  "bad"  tape. 
Replace  this  tape  with  another  reel,  and  press  START  key.  If  this  STOP 
persists,  transfer  via  console  to  location  (213)g  to  by-pass  the  dump.  The 
program  will  not  be  affected  unless  a  subsequent  machine  error  occurs,  in 
which  case  it  will  be  necessary  to  restart  the  program  from  the  beginning. 

With  the  programs  stored  permanently  on  magnetic  tape,  the  data  cards 
will  be  read  via  the  on-line  card  reader  as  follows; 

1.  Load  the  FORTRAN  START  card. 

2.  After  IB  Monitor  begins  to  function  (a  few  seconds),  FEED  out  this 
START  card. 

3.  Place  data  cards  in  card  reader  hopper  and  depress  reader  START 
key  until  READY  light  goes  on. 

4.  Program  will  automatically  read  these  cards,  for  END-OF-FILE 
convenience,  place  two  blank  cards  at  the  end  of  the  deck  of  data  cards. 
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VI.  COMPARISON  OF  RESULTS 

.  -  Vibration  Program: 

The  table  below  lists  the  experimental  frequencies  obtained  from  a  model 
of  the  DD  828  propeller  blade  (Reference  4),  and  computed  frequencies  utilizing 
45  degrees  of  freedom,  which  was  the  limit  for  the  previous  IBM  704  program, 
and  90  degrees  of  freedom.  An  8x9  quadrature  integrating  mesh  was  used  for 
both  computer  runs.  A  density  of  0.  31  lbs.  /cu.  in.  and  a  modulus  of  elasticity 
of  15  •  10^  psi,  were  assumed.  The  scale  factor  between  the  model  ?.nd  the  actual 
blade  is  10. 125. 

The  comparison  is  made  on  the  basis  of  similar  mode  shapes.  Note  that 
the  mode  shap^  for  the  third  experimental  mode  corresponds  with  that  obtained  by 
the  computer  for  its  second  mode,  and  vise-versa. 

MODEL 


Mode  Computed  Experimental 


45 

90 

1 

573 

480 

446 

2 

949  ' 

800 

981  (3) 

3 

1140 

1017 

945  (2) 

4 

1392 

1139 

1220 

5 

1496 

1392 

-- 

6  ■ 

2005 

1'610 

1425 

7 

2441 

1845 

1640 

8 

2577 

2111 

1890 

9 

2753 

2210  (84) 

2097 

10 

3062 

2450  (84) 

2350 
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The  tables  on  the  following  pages  indicate  the  asymptotic  behavior  of  the 
eigenvalues  (and  natural  frequencies),  as  the  number  of  degrees  of  freedom 
prescribed,  is  increased.  This  behavior  is  graphically  illustrated  on  page  43. 

The  mode  shapes  are  shown  on  pages  44  to  53  ,  The  solid  lines  are  drawn 
from  data  generated  by  the  program;  the  broken  lines  were  those  detected 
experimentally.  It  became  increasingly  difficult  for  the  experimenter  to  obtain 
the  mode  shapes  as  the  frequency  increased;  therefore,  no  comparison  is 
attempted  beyond  mode  number  7. 


COMPARISON  OF  EIGENVALUES 


U1 

CO 

o 

o 

o 

o 

o 

o 

o 

O 

p^ 

p-^ 

pH 

pH 

• 

• 

• 

• 

• 

• 

•. 

o 

00 

00 

o 

03 

00 

Tti 

rjt 

r- 

00 

sO 

fs3 

4* 

(S) 

sD 

Oi 

CO 

in 

r- 

ON 

• 

• 

• 

• 

• 

• 

• 

• 

O 

o 

o 

o 

o 

o 

o 

o 

m 

T}* 

CO 

(O 

CO 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

^H 

ii 

r» 

(O 

d 

Co 

d 

00 

00 

o* 

00 

m 

CO 

CO 

sD 

sD 

CO 

CO 

m 

0- 

iH 

O 

o 

o 

o 

o 

o 

o 

d 

d 

d 

iri 

Tt* 

■tj' 

CO 

CO 

CO 

CO 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

m 

rH 

rH 

• 

• 

• 

• 

• 

• 

• 

• 

• 

CO 

'X) 

p^ 

sD 

00 

CO 

o 

CO 

O' 

r* 

o 

o 

00 

o 

CO 

CO 

r* 

CO 

CO 

m 

0- 

pH 

pH 

pH 

rH 

o 

o 

o 

o 

o 

o 

d 

d 

o 

d 

LO 

co 

CO 

CO 

CO 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

pH 

fH 

pH 

rH 

ro 

-h 

00 

in 

b^ 

■4* 

CO 

r- 

o 

04 

in 

o 

r- 

<^J 

CO 

-t> 

in 

00 

»H 

pH 

rsj 

o 

o 

o 

o 

o 

o 

d 

d 

d 

o 

m 

4* 

-><< 

fO 

CO 

CO 

CO 

CO 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

p^ 

•H 

pH 

rH 

rH 

in 

<M 

o 

ro 

P^3 

vO 

in 

pH 

pH 

m 

nO 

rsj 

^o 

00 

rH 

o 

O' 

rvj 

CO 

in 

sD 

i-H 

pH 

n3 

(M 

o 

o 

o 

o 

O 

d 

d 

o 

o 

O 

•>t 

■>i< 

fO 

fO 

CO 

CO 

CO 

o 

o 

o 

o 

O 

o 

o 

o 

O 

O 

f-H 

pH 

1^ 

rH 

rH 

pH 

o 

m 

00 

in 

in 

in 

pH 

sC 

in 

rsi 

m 

r- 

r- 

pH 

CO 

o 

CO 

CJv 

o 

f*4 

^3 

CO 

o 

o 

o 

o 

o 

d 

O 

d 

o 

d 

41 


44 


mode  no.  I 


mode  no.  l 


propeller  blade 


Page 


4- 


PROP 
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Stress  Program 

The  stress  program  is  subject  to  the  same  round-off  error  problems 
associated  with  the  vibration  program.  Indeed,  the  triple-precision  technique 
used  to  factor  the  strain  energy  matrix,  is  identical  for  the  two  programs. 
Thus  the  same  restrictions  on  maximum  size  of  the  system  that  apply  for  the 
vibration  program,  also  apply  for  the  stress  program, 

A  parallel  series  of  test  runs  was  conducted  on  the  DD828  propeller 
blade,  using  a  "peaked"  pressure-jump  loading.  Two  runs  were  conducted 
at  84  degrees  of  freedom,  one  using  twice  as  many  iterations  to  improve 
the  solution  vector,  as  is  normally  used.  The  results  of  the  two  runs 
were  practically  identical,  indicating  that  an  adequate  number  of  iterations 
are  being  used  and  that  the  solution  of  the  system  of  linear  equations,  is 
accurate.  Another  series  of  test  runs  was  conducted  on  a  stainless  steel 
propeller,  loaded  with  a  uniform  pressure  distribution. 

Discussion 

The  results  of  the  first  series  of  runs  were  conducted  using  45,  63  and 
84  degrees  of  freedom,  and  are  indicated  graphically  on  the  following  pages. 

In  addition,  detailed  displacement  and  stress  distributions  are  plotted  over 
the  blade  planform  for  the  run  utilizing  84  degrees  of  freedom. 

It  is  worth  noting  that  as  the  size  of  the  system  grows,  the  free  boundary 
conditions  are  satisfied  with  increasing  accuracy.  The  graph  on  page  57 
clearly  indicates  this  trend,  with  the  span-wise  stress  (d**)  virtually  vanishing 
at  the  blade  tip  for  the  run  utilizing  84  degrees  of  freedom.  In  addition,  this 
trend  is  also  true  for  the  chord-wise  stress  ((T**V  indicated  on  page  56, 

The  final  results  obtained  by  using  84  degrees  of  freedom  for  the  stainless 
steel  propeller,  are  graphically  shown  on  pages  58  to  67  , 
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vn  CONCLUSIONS 

An  exhaustive  series  of  computer  test  runs  were  conducted  on  a  particular 
propeller  blade  configuration  (DD  828),  which  was  considered  to  be  representa¬ 
tive  of  blades  characterized  by  severe  built-in  twist,  extremely  low  aspect 
ratio,  and  large  overhang  of  thin  leading  and  trailing  edges.  The  results  of 
the  runs  made  by  the  vibration  program  indicate  conclusively  that  it  is  necessary 
to  prescribe  a  large  number  of  degrees  of  freedom  to  obtain  satisfactory  results, 
for  this  type  of  blade.  Unfortunately,  the  marked  growth  of  accumulated  round¬ 
off  error  as  the  size  of  the  system  increases,  presents  a  serious  problem. 
Indeed,  this  factor  defines  the  upper  bound  of  the  allowable  number  of  degrees 
of  freedom  which  may  be  prescribed. 

This  IBM  7090  program  was  designed  to  accommodate  up  to  108  degrees 
of  freedom,  whuch  corresponds  to  a  production  machine  time  of  approximately 
1-1/Ehr.s.This  limit  was  imposed  by  the  storage  capacity  of  the  computer,  as 
well  as  by  economic  considerations  with  respect  to  cost  of  machine  time 
necessary  to  complete  a  single  run.  It  was  believed  that  this  number  would 
be  adequate  to  cause  the  natural  frequencies  to  "settle  down"  and  behave 
asymptotically. 

The  accumulated  round-off  error  did  not  permit  this  series  of  runs  to 
attain  the  maximum  size  for  which  the  program  is  designed.  It  was  necessary 
to  resort  to  triple-precision  arithmetic  in  a  part  of  the  eigenvalue  solution,  to 
'  curb  the  rovind-off  error  sufficiently  to  allow  a  successful  run  at  90  degrees 

of  freedom.  Despite  this  refinement,  the  accumulated  round-off  error 

% 

caused  extraneous  negative  values  for  2  eigenvalues  to  be  printed  in  the  output. 
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These  extraneous  values  apparently  did  not  effect  the  remaining  88  eigenvalues 
as  examination  of  pages  to  will  verify;  they  did  indicate  that  the  accumulated 
errors  were  becoming  critical,  and,  as  a  matter  of  fact,  a  subsequent  trial 
run  using  96  degrees  of  freedom  failed. 

It  is  believed  that  the  bulk  of  the  accumulated  round-off  errors  is  generated 
in  the  process  of  forming  and  integrating  the  energy  matrices.  To  alleviate 
diis  condition,  it  would  be  necessary  to  utilize  double -precision  arithmetic  iii 
this  section  of  the  program.  This,  however,  would  be  a  major  programming 
effort  and  the  increased  cost  of  machine  time  could  triple. 

An  8x9  integrating  quadrature  mesh  is  used.  Test  runs  made  with  10  x  10 
and  12x12  meshes  exhibited  no  significant  changes  in  the  eigenvalues;  thus, 
these  expensive  refinements  do  not  contribute  to  any  increase  in  accuracy. 

For  the  great  majority  of  blade  configurations,  it  is  believed  that  the 
present  program  will  generate  accurate  estimates  of  the  natural  frequencies. 

This  series  of  runs  indicates  that  the  natural  frequencies  have  indeed  "settled 
down"  at  90  degrees  of  freedom,  and  that  these  figures  can  be  referred  to, 
with  confidence.  The  mode  shapes  remained  essentially  unchanged  between 
75  and  90  degrees  of  freedom. 

The  limitation  on  size  due  to  the  accumulation  of  round-off  error,  represents 
a  much  more  severe  restriction  on  stress  program  than  on  the  vibration  program. 
The  natural  modes  of  vibration  are  an  "integrated"  property  of  the  entire  blade 
and  thus  are  intrinsically  more  responsive  to  "settling  down"  as  the  size  of  the 
system  grows.  On  the  other  hand,  the  stress  program  mesh  furnish  quantitatively 
precise  data  over  the  entire  domain  of  the  blade  planform.  To  attain  the  same 
degree  of  accuracy  for  the  stress  distributions  as  is  noW  obtained  for  the 


vibration  frequencies,  if  is  necessary  to  add  more  Rayleigh  functions, 
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which  leads  to  larger  systems  of  equations. 

The  method  now  used  to  solve  the  system  of  linear  equations,  depends 
on  this  system  being  positive -definite.  This  condition  is  negated  by  the  round¬ 
off  error  incurred  in  forming  a  large -size  strain  energy  matrix;  a  technique 
to  solve  the  system  of  equations  which  does  not  depend  upon  the  property 
of  positive  definiteness  must  be  employed  in  order  to  accommodate  a 
size  larger  than  say  84  by  84  .  A  Gaussian  elimination  process,  followed 
by  a  Seidelian  iteration  to  improve  the  solution  vector  should  accomplish 
this  result. 

The  cost  of  running  the  present  stress  program  has  been  reduced  from 
1-1/2  hours  to  approximately  55  minutes.  See  Appendix  II  for  details. 

This  cost  could  be  further  reduced  25-30%  by  converting  the  equipment  to 
an  IBM  7094  with  high-speed,  high-density  tape  units. 
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The  Rayleigh  functions  must  be  chosen  to  satisfy  the  displacement  boundary 

aaj 

conditions  at  the  root,  vis.  u  =  u  =  0  and  -  =  u  =  0,  The  geometry  of  the 

*  *  1  ^ 

root  is  approximated  by  a  parabola  of  the  form,  y®^*ay^^  +  by^*^+  c.  Since 

the  root  passes  through  the  origin,  the  constant  term,  c,  vanishes.  The 
coefficients  a,  and  b,  are  expressed  in  terms  of  the  coordinates  of  the  end¬ 
points  of  the  root,  which  are  given  on  input  card  2. 

The  boundary  conditions  for  u  and  u  are  satisfied  by  the  functions,  , 
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APPENDIX  II 

A  great  variety  of  techniques  have  been  developed  and  enripioyed  to 
solve  large  systems  of  simultaneous  equations.  The  relative  usefulness  of 
each  technique  is  intimately  related  to  the  size  and  character  of  the  system  in 
question..  At  all  times,  the  programmer  must  contend  with  the  problem  of 
propagated  round-off  errors  and  recognize  the  economic  necessity  of 
minimizing  machine  time,  consistant  with  the  desired  accuracy. 

In  the  course  of  solving  the  two -matrix  eigenvalue  problem  it  was 
necessary  to  factor  the  strain  energy  matrix  using  multiple -precision 
arithmetic  which  utilized  two  "words"  for  the  fraction  -  a  total  of  71  bits 
compared  with  25  bits  used  for  single  precision  (plus  sign).  After  attemping 
several  other  techniques  to  solve  the  system  of  equations  in  the  stress 
program  all  of  which  proved  to  be  unacceptable  for  one  reason  or  another,  it 
was  decided  to  utilize  the  existing  coding  of  the  factoring  technique  to  generate 
an  approximate  inverse,  and  then  try  to  improve  the  solution  by  various 
means.  Three  techniques  are  used,  in  series;  the  improvement  effected 
by  one  procedure,  materially  assists  the  convergence  of  the  following  one. 

1.  Series  Expansion  to  Improve  the  Inverse 

Given  a  matrix.  A,  and  an  approximate  inverse  A  ^=A  ^  +E,  then 
A  a"^=  A(A'^+E)  =  I  +  C  where  C  =  AE 
The  matrix  C  may  be  calculated;  E  is  the  error  in  the  approximate  inverse  . 
Solving  for  E,  E  =  A'^C 

E  -f  EC  =  E(I+C)  =  a'^C+EC  =  (A'^+E)C  =  A  ’^C 
E(I+C)(I-C)  =E(I-C*)  =  A'^C(I-C) 

Multiplying  by  (I+C*);  EU-C^)  =  a"\c-C*+C*-C^) 
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In  general,  E(I-C®^)  =  a"^(C+C*+C^  - _ -  c”) 

Thus  the  error  matrix,  E,  may  be  calculated,  subject  to  a  truncation  error, 
-EC^.  If  the  matrix  C  is  "sufficiently”  small,  the  series  expansion  will 
converge,  and  the  truncation  error  will  be  extremely  small,  compared  to  E. 
Then,  since  A”^  =  A“^  -  E,  we  may  write 

a’^  =  A"  -  a  (C-C^+C^  -  . . ,  -  c”)  or 

A-1  =  A  (I-C+C^-C^  ...  +  C’^) 

For  this  program,  n  =  8,  If  C  is  "too  large",  the  series  will  diverge 
amd  the  entire  improvement  of  the  solution  will  be  done  by  the  Seidclian 
iteration  technique.  Following  this  series  expansion(if  successful)  is: 

2.  Iteration  on  Forcing  Vector 

Given  the  matrix  A,  a  good  approximation  to  its  inverse  ,  A“^,  and  the 

forcing  vector,  p,  in  the  system  Aq  =  p.  Then  we  may  write 
~-l 

qj^  =  A  p  where  q^^  is  the  first  approximation  to  the  solution  vector. 

To  improve  it,  compute  p^;  Aq^  =  p^^  or  A(q  +  Aq^)  =  p+Ap^  .  Then  we 

may  solve  for  the  error  in  the  solution  vector,  Aqj^  =  A  ^Apj^  and  writing 

qj  =  qi  +  A  qj,  can  repeat  the  process.  For  this  iteration  to.  converge,  the 
~  -1 

inverse,  A  ,  must  be  "sufficiently"  accurate.  Otherwise,  we  must  depend 
entirely  upon  the  Seidelian  iteration, 

3.  Seidel  Iteration  on  the  Solution  Vector 

Transposing  the  i^^  equation  of  the  system,  we  may  solve  for  the  i^^ 
element  of  the  solution  vector  in  terms  of  the  best  values  computed  so  far. 


for  the  other  n  -  1  elements  of  the  solution  vector: 
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[pi  -  C  ‘ij  ’j] 

j=Ai 


The  program  sweeps  through  the  system  many  times  until  either  a 
satisfactory  tolerance  is  achieved  or  until  the  upper  limit  of  iterations 
is  attained.  Since  this  iteration  technique  is  unconditionally  convergent 
for  positive  definite  symmetric  matrices,  the  final  solution  vector  will 
have  the  necessary  accuracy.  The  upper  limit  of  the  number  of  iterations 
depends  upon  the  success  attained  by  the  previous  two  techniques.  This 
upper  limit  is  significantly  reduced  (thus  saving  machine  time)  if  the 
prior  iterations  were  successful  to  some  extent. 
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Fig.  1 

Blade  Geometry  and  Coordinate  Syetems 
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